i = 0
gbm.results <- list()
ntree_seq <- seq(5,150,5)
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Fit GBM
  gbm.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("gbm",
                                       "pROC")) %dopar% {   # Loop over CV runs
                           set.seed(52*cv) 
                           shuffle <- sample(1:N,
                                             N,
                                             replace=FALSE)
                           pred.prob <- matrix(NA,nrow=N,ncol=length(ntree_seq))
                           sampleSplit = c(0,round((N/5*c(1:5))))
                           for (f in 1:5) {
                             # Fit lasso
                             cv.fit = gbm.fit(x = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                        rhs]),
                                              y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                          dv]),
                                              n.trees = max(ntree_seq),
                                              distribution = gbm_params$distribution,
                                              interaction.depth = gbm_params$interaction.depth,
                                              n.minobsinnode =gbm_params$n.minobsinnode,
                                              verbose = gbm_params$verbose,
                                              shrinkage = gbm_params$shrinkage,
                                              bag.fraction=gbm_params$bag.fraction,
                                              train.fraction=gbm_params$train.fraction
                             )
                             pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                       ] <- predict(cv.fit,
                                                    newdata = as.matrix(dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                      rhs]),
                                                    type = "response",
                                                    n.trees = ntree_seq)
                           }
                           cv_results <- c("CV Run"=0,
                                           "Trees"=0,
                                           "MSE"=1000)
                           # Get error rates for each lambda used
                           for (t in 1:length(ntree_seq)) {
                             if (is.na(pred.prob[1,t])==0) {
                               ###!####
                               mse <- mean((dta.past[,dv] - pred.prob[,t])^2)
                               ### ! ####
                               if (mse<cv_results["MSE"]) {
                                 cv_results <- c("CV Run"=cv,
                                                 "Trees" = ntree_seq[t],
                                                 "MSE" = mse)
                                 best.t <- t
                               }
                             }
                           } 
                           cv_results<- c(cv_results,
                                          "DT" = as.numeric(mostaccDT(y=dta.past[,dv],
                                                                      yhat=pred.prob[,best.t],
                                                                      J=1000)),
                                          "DTsens" = as.numeric(bestDT(y=dta.past[,dv],
                                                                      yhat=pred.prob[,best.t],
                                                                      J=1000)))
                           return(cv_results)
                  }
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  if (cv.runs>1) {
    gbm.results[[i]]$trees <- round(mean(cv_results[,"Trees"]))
    gbm.results[[i]]$dt <- mean(cv_results[,"DT"])
    gbm.results[[i]]$dtsens <- mean(cv_results[,"DTsens"])
  } else {
    gbm.results[[i]]$trees <- cv_results["Trees"]
    gbm.results[[i]]$dt <- cv_results["DT"]
    gbm.results[[i]]$dtsens <- cv_results["DTsens"]
  }
  set.seed(i)
  mod = gbm.fit(x = as.matrix(dta[dta[,t.var]< year,
                                  rhs]),
                y = as.matrix(dta[dta[,t.var]< year,
                                  dv]),
                n.trees = gbm.results[[i]]$trees,
                distribution = gbm_params$distribution,
                interaction.depth = gbm_params$interaction.depth,
                n.minobsinnode =gbm_params$n.minobsinnode,
                verbose = gbm_params$verbose,
                shrinkage = gbm_params$shrinkage,
                bag.fraction=gbm_params$bag.fraction,
                train.fraction=gbm_params$train.fraction
  )
  gbm.results[[i]]$fit.oos <- predict(mod,
                                      newdata = as.matrix(dta[dta[,t.var] == year,
                                                                   rhs]),
                                      type = "response",
                                      n.trees = gbm.results[[i]]$trees)
  gbm.results[[i]]$actual.oos <- as.matrix(dta[dta[,t.var] == year,
                                                    dv])
  gbm.results[[i]]$param_dist <- cv_results
  
  print(year)
}
save(gbm.results,file=paste(modeldir,"/",
                              country,
                              "_gbm_",
                              v,
                              "_",
                              rhs.group,
                              "_mse",
                              ".RData",
                              sep = ""))